General Practitioners and population by Medicare Primary Health Network area

### start necessary libraries

# data import libraries
library(xlsx)
library(foreign)
library(rgdal)

library(tidyverse)
library(dplyr)

# reproducible research library
library(knitr)

# plotting and maps libraries
library(rmapshaper) # for simplifying maps
library(ggplot2)
library(leaflet)
library(ggmap)

### set figure save path to figure/
knitr::opts_chunk$set(fig.path = 'figure/')
knitr::opts_chunk$set(fig.width=12, fig.height=8)

Download and import data files

Population data in Australia 2009 to 2014, grouped by Statistical Area Level 3 (SA3) from http://www.health.gov.au/internet/main/publishing.nsf/Content/PHN-Demographic_Data

Primary Health Network (‘PHN’) map data from http://www.health.gov.au/internet/main/publishing.nsf/content/phn-boundaries

Concordance table between PHN geographical area and Statistical Area Level 3 from http://www.health.gov.au/internet/main/publishing.nsf/Content/PHN-Concordances

General Practice (‘GP’, also known as family medicine) workforce by PHN found from Australian Government Health Workforce Data website (http://hwd.health.gov.au/), direct links for General Practitioner data by PHN is not available without login and setting search terms. Pre-downloaded 2016 data.

## OGR data source with driver: ESRI Shapefile 
## Source: "PHN_boundaries_AUS_May2017_V7.shp", layer: "PHN_boundaries_AUS_May2017_V7"
## with 31 features
## It has 8 fields
## Integer64 fields read as strings:  OBJECTID

Place information into Map file

Places Population, General Practitioner numbers and Population per General practitioner information into map file

Creates a data frame with

# choose all ages statistic
select.age.population <- select(population.2014, SA3_NAME_2011, AllAges)

# calculate population in each PHN
phn.population <- sa3.phn.concordance %>%
  merge(select.age.population, by = 'SA3_NAME_2011') %>%
  # merge population numbers by SA3 areas
  mutate(PHN_Number = as.numeric(AllAges)*as.numeric(RATIO)) %>%
  # calculate population in each SA3 as proportion in each PHN
  # some SA3 are spread across more than one PHN
  group_by(PHN_NAME) %>%
  # PHN can have multiple SA3
  summarize(population = sum(PHN_Number)) 

# add GP information to each PHN
# calculate population per GP ratio in each PHN
combined.oz.phn.description <- oz.phn.description %>%
  merge(gp.numbers, by = 'PHN_NAME') %>%
  merge(phn.population, by = 'PHN_NAME') %>%
  mutate(pop.gp.ratio = population/gp) %>%
  # just select the variables that we need
  select(PHN_NAME, gp, population, pop.gp.ratio)

add extra map data to the map

oz.phn.map@data <-  merge(oz.phn.map@data, combined.oz.phn.description, by = 'PHN_NAME')

Population and General Practitioner numbers in Australia, by Primary Health Network areas

‘Click’ on region to reveal Population and GP numbers

oz_loc <- geocode('Australia', source='dsk')
## Information from URL : http://www.datasciencetoolkit.org/maps/api/geocode/json?address=Australia&sensor=false
phn.map <- ms_simplify(oz.phn.map) %>% # reduce points to about 5% of original
  leaflet() %>% # the map generator
  addTiles() %>% # base map, usually OpenStreetMap
  setView(lng = oz_loc$lon, lat = oz_loc$lat, zoom = 4) %>% # set view
  addPolygons(color="#444444", weight=1, smoothFactor=0.5,
              opacity=0.8, fillOpacity = 0.15,
              popup = ~paste(PHN_NAME, '<br>Population :',
                             sprintf('%.0f',population),
                             '<br>GPs :', gp,
                             '<br>Population to GP :',
                             sprintf('%.0f',pop.gp.ratio)),
              fillColor = ~colorQuantile("YlOrRd", pop.gp.ratio)(pop.gp.ratio),
              highlightOptions =
                highlightOptions(color='green',weight=2,bringToFront = TRUE))

phn.map

‘Click’ on region to reveal Population and GP numbers